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Abstract 

We generalize the Euclidean 1-center approximation algorithm of [Badoiu and 



Clarkson (2003) to arbitrary Riemannian geometries, and study the corre- 
sponding convergence rate. We then show how to instantiate this generic 
algorithm to two particular settings: (1) the hyperbolic geometry, and (2) 
the Riemannian manifold of symmetric positive definite matrices. 

Keywords: 1-center; minimax center; Riemannian geometry; core-set; 
approximation 



1. Introduction and prior work 



Finding the unique smallest enclosing ball (SEB) of a finite Euclidean 
point set P = {pi,...,p„} is a fundamental problem that was first posed 
by Sylvester ( |1857 ). This problem has been thoroughly investigated in the 



computational geometry community by Welzl (1991) and Nielsen and Nock 



(2009), where it is also known as the minimum enclosing ball (MEB), the 
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1-center problem, or the minimax optimization problem in operations re- 
search. In practice, since computing the SEB exactly is intractable in high 
dimensions, efficient approximation algorithms have been proposed. An al- 



gorithmic breakthrough was achieved by Badoiu and Clarkson (2008) that 
proved the existence of a core- set C C P of optimal size \C\ = P]To that 
< (1 + e)r{P) (for any arbitrary e > 0), where r{S) denotes the radius 
of the SEB of S. Let c{S) denote the ball center, i.e. the minimax center. 
Since the size of the core-set depends only on the approximation precision e 
and is independent of the dimension, core-sets have become widely popular 
in high-dimensional applications such as supervised classification in machine 



learning (see for example, the core vector machines of Tsang et al. (2007)). 



In the work of Badoiu and Clarkson (2003), a fast and simple approximation 



algorithm is designed as follows: 



BC-ALG: 



Initialize the center ci G P, and 

Iteratively update the current center using the rule 

fi Cj 



Ci+i ^ Ci + 



i + l 



where /j denotes the farthest point of P to Cj 



It can be proved that a (l + e)-approximation of the SEB is obtained after 
[^] iterations, thereby showing the existence of a core-set C = {/i,/2,...} 
of a size at most [^]: r(C) < (1 + e)r{P). This simple algorith m runs in 
time O(^), and has been generalized to Bregman divergences by 



Nock and 



Nielsen (2005) which include the (squared) Euclidean distance, and are the 
canonical distances of dually flat spaces, including the particular case of self- 
dual Euclidean geometry. (Note that if we start from the optimal center ci = 
c{S), the ffist iteration yields a center C2 away from c{S) but it will converge 



in the long run to c{S).) Badoiu and Clarkson (2008) proved the existence 
of optimal e-core-set of size [-]. Since finding tight core-sets requires as a 
black box primitive the computation of the exact smallest enclosing balls of 
small-size point sets, we rather consider the Riemmanian generalization of 
the BC-ALG, although that even in the Euclidean case it does not deliver 
optimal size core-sets. 



2 



Many data-sets arising in medical imaging (see Pennec (2008)) or in com 



puter vision (refer to Turaga and Chellappa (2010)) cannot be considered as 



emanating from vectorial spaces but rather as lying on curved manifolds. For 
example, the space of rotations or the space of invertible matrices are not 
fiat, as the arithmetic average of two elements does not necessarily lie inside 
the space. 

In this work, we extend the Euclidean BC-ALG algorithm to Riemannian 
geometry. In the remainder, we assume the reader familiar with basic notions 



of Riemannian geometry (see Berger (2003) for an introductory textbook) in 
order not to burden the paper with technical Riemannian definitions. How- 
ever in the appendix, we recall some specific notions which play a key role in 
the paper, such as geodesies, sectional curvature, injectivity radius, Alexan- 
drov and Toponogov theorems, and cosine laws for triangles. Furthermore, 
we consider probability measures instead of finite point set^ so as to study 
the most general setting. 

Let M be a complete Riemannian manifold and v a probability measure 
on M. Denote by p(x, y) the Riemannian distance from x to ?/ on M that 
satisfies the metric axioms. Assume the measure support supp(z/) is included 
in a geodesic ball Bio, R). 

Recall that if p G [1, oo) and / : M — )■ M is a measurable function then 



LPiu) 



M 



1/p 



and 



Let 



inf{a>0, u{{yeM, \f{y)\>a}) = 0} 



■a,p 



1 



mm 
min 



in{inj(M),^} if l<p<2, 
{inj(M),^} if 2<p<oo 



(1) 



where inj(M) is the injectivity radius (see the appendix) and a > is such 
that is an upper bound for the sectional curvatures in M (in fact replacing 
M by B{o, 2R) is sufficient, so that we can always assume that a > 0). For 
p G [1, C)o], under the assumption that 



R<R. 



a,p 



(2) 



^We view finite point sets as discrete uniform probability measures. 
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it has been proved by Afsari (2011) that there exists a unique point Cp which 
minimizes the following cost function 

Hp : M [0, oo] 

X ^ \\p{x,-)\\lp{u) (3) 

with Cp G B{o, R) (in fact, lying inside the closure of the convex hull of the 
support of z/). 

For a discrete uniform measure viewed as a "point cloud" in an Euclidean 
space and p G [l,C)o), we have Hp{x) = {^Yl^=i \\Pi ~^llp)^^^5 with || • ||p 
denoting the Lp norm, and Hoc (x) is the distance from x to its farthest point 
in the cloud. 

In the general situation the point Cp that realizes the minimum represents 
a notion of centrality of the measure (eg., median for p = 1, mean for p = 2, 
and minimax center for p = oo). This center is a global minimizer (not only 
in B{o,R)), and this explains why a bound for the sectional curvature is 



required on the whole manifold M (in fact B{o,2R) is sufficient, see Afsari 
( |2UIT| )). 

Deterministic subgradient algorithms for finding Cp have been considered 



by Yang (2010) for the median case {p = 1). Stochastic algorithms have been 



investigated by Arnaudon et al. ( 2010[ ) for the case p G [1, oo), and a central 



limit theorem (CLT) for the suitably renormalized process is derived (in fact 
a convergence in law to a diffusion process). See also for similar algorithms 



minimizing other cost functions, the work of Bonnabel (2011). 

In this work, we consider the case p = oo, with Coo denoting the minimax 
center. Hereafter we use c for Coo, H for H^o and Ra for Ra,oo- In this case 
there is no canonical deterministic algorithm which generalizes the gradient 
descent algorithms considered for p G [1, oo). Following Eq. [sj H{x) denotes 
the farthest distance from x to a point of the support of the measure (L°°- 
norm) . 

To give an example of a Riemannian manifold, consider the space of 
symmetric positive definite matrices with associated Riemannian distance 
(see Section |4]) 



p{p,Q) = ||iog(p-ig)||^ = (4) 

where Aj are the eigenvalues of matrix P~^Q. This is a non-compact Rieman- 
nian symmetric space of nonpositive curvature (Cartan-Hadamard manifold. 



4 



see Lang (1999), chapter 12). In this context any measure u with bounded 
support satisfies. Eq. [2] (since we can take a > as small as we like), and con- 
sequently the minimizer c of H exists and is unique. We call it the 1-center 
or minimax center of u. 

We generalized the BC-ALG by noticing that the iterative update is a 
barycenter of the current minimax center with the current farthest point. 
Thus the new position of the minimax center falls along the straight line 
joining these two points in Euclidean geometry. In Riemannian geometry, 
the shortest path linking two points is called a geodesic (for example, arc of 
a great circle for spherical geometry). Instead of walking on a straight line, 
we instead walk on the geodesic to the farthest point as follows: 



GEO-ALG: 



Initialize the center with ci G P, and 
Iteratively update the current minimax center as 

1 



Q+i = Geodesic q, fi, . , ^ 
V « + 1 

where fi denotes the farthest point of P to q, and 
Geodesic(j9, q, t) denotes the intermediate point m on the 
geodesic passing through p and q such that p(j9, m) = tx p{p, q). 



Note that GEO-ALG generalized BC-ALG by taking the Euclidean dis- 
tance p{p, q) = \\p — q\\. 

The paper is organized as follows: Section [2] gives and proves a crucial 
lemma. It is followed by the description and convergence rate analysis of 
our generic Riemannian algorithm in Section |3} Section |4] instantiates the 
algorithm for the particular cases of the hyperbolic manifold and the manifold 
of symmetric positive definite matrices. Section |5] concludes the paper and 
hints at further perspectives. To make the paper self-contained, the appendix 
recalls the fundamental notions of Riemannian geometry used throughout the 
paper. 
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2. A key lemma 

In this section, we assum^that supp(z/) C B{o,R) and 

1 . n 



R<R^ = -min{inj(M),-| 
2 L a J 



with a > such that is an upper bound for the sectional curvatures in M. 
The following lemma is essential for proving the convergence of the algorithm 
determining the minimax of u. 

Lemma 1. There exists r > such that for all x G B{o,R), 

H{x) - H{c) >Tp'^{x,c). (5) 

Proof: 

The point c is the center of the smallest ball which contains supp(z/) and the 
radius of this ball is exactly r* := H{c) (see [Afsarl (2009)). An immediate 



consequence is that r* < R. Denoting by S{c,r*) the boundary of this ball 
and by ScM the set of unitary vectors in T^M, for all v G ScM there exists 
y G S{c,r*) n supp(z/) such that 

{c^,v)<0 (6) 

where t i— )■ 'yt{c,y) is the geodesic from c to y in time one, %{c,y) denotes 
derivative with respect to t and = 7o(c, y). Indeed, if this was not true it 
would contradict the minimality of S'(c, r*) (refer to Afsari ( 2009[ )). 



Now letting t 7j(f ) = exp2.(tf) the geodesic satisfying 70(f) = f, we 
prove Eq. [slfor x = 7t(f ). We have 



Hiltiv)) -H{c) > p{-ft{v),y) - p{c,y) = p{-ft{v),y) - r* (7) 

by definition of H. 

Then we consider a 2-dimensional sphere with constant curvature 
a^, distance function p, and in 5*^2 a comparison triangle '^t{v)yc such that 
pijj, c) = r*, f is a unitary vector in TgS'^a satisfying 

dy,v) = {c^,v) (8) 



^Any bounded measure on a Cartan-Hadamard manifold satisfies this assumption. 



6 



Let us prove that 

p{lt{v),y) -r* = p{^t{v), y) - p{c, y) > T^p^{^t{v), c) (9) 

for some Tq, > provided condition Eq. |6] is realized: for simplicity we will 
write d = p{^t{v),y). Using Eq. ojand the first law of cosines (Theorem [i] in 
the appendix), we get 



\ cos f arf j — cos (ar*) cos(at) 

cos I cy ^ V \ = ^ — 

~ V ' / sin (ar*) sin(Q;t) 

which yields 

cos {ad^ — cos (ar*) cos(at) < 0. 
On the other hand since < ar* < | and < ad < tt we have 



(10) 



< 2 sin \ adj cos(ar*) sin(Q;r*). 

So we get 

cos ^adj — cos (ar*) cos(at) < 2 sin {ad^ cos(ar*) sin(ar*) 

which is equivalent to 

cos^(ar*) cos(a(i) + cos(ar*) sin(ar*) sin(a(i) — cos(ar*) cos(at) 
< sin(ac/) cos(ar*) sin(ar*) — sin^(ar*) cos(ac/) 

and this in turn implies 

sin ^a — r* j j > cotan(ar*) ^cos ^a — r* j j — cos(at)^ 

so 

. rP{ltiv),y)-r* a *n ^ « . / 

hmmt > — cotan(ar ) > — cotan(aitaj 



uniformly in v. Consequently Eq. |9]is true for 7t('u) in a neighborhood of c, 
and since p{^i{v)^y) — r* does not vanish outside this neighbourhood, by a 
compactness argument we prove that Eq. |9] is true in any compact included 
in B{c, Ra), if Tq, is sufficiently small. 

To finish the proof we are left to use the Alexandrov comparison theorem 
(Theorem |2]in the appendix) with triangles 'yt{v)yc and jtiv)yc to check that 
the right hand side of Eq. [7] in M is larger than the left hand side of Eq. |9] 
This proves Eq. [s] in B(c,R) fl B(o,R), and for proving it in B{o,R) we 
just have to notice that H is continuous and positive on the compact set 
B{o, R)\B{c, R), hence it has a positive lower bound. □ 
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3. Riemannian approximation algorithm 

For X G B{o,R), denote by t i— 7t(t>(x, z/)) a unit speed geodesic from 
7o(f (x, z/)) = X to one point y = 'yH{x){v{x, z/)) in supp(z/) wliicli realizes the 

maximum of the distance from x to suppfz/). So f = , exp~^(y). A 

H{x) 

measurable choice is always possible. Note that if v has finite support, when 
there is a finite number of possibilities for y it is natural to make a random 
uniform choice. However in a generic situation this should never happen, 
there should be only one choice. 

We consider the following stochastic algorithm. 



RIE-ALG: 
Fix some 5 > 0. 

Step 1 Choose a starting point Xq G supp(z/) and let /c = 

Step 2 Choose a step size t^+i ^ (0, 5] and let Xk+i = Itk+ii'^i^k, ^)), 

then do again step 2 with A; A; + 1. 



This algorithm generalizes the Euclidean scheme of Badoiu and Clarkson 



(2003) and algorithm GEO-ALG for probability measures. Indeed, if GEO- 
ALG is initialized with G P with the first integer larger than 1/6, then 
it suffices to take tk = l/kioi k> ko in RIE-ALG. 

Let a Ab denote the minimum operator a Ab = min(a, b). 

Let 

iJo = ^A-. (11) 

Theorem 1. Assume a,P > are such that — is a lower bound and o? 
an upper bound of the sectional curvatures in M. 
If the step sizes (tk)k>i satisfy 

R 2 

S < ^ A -arctanh (tanh(^i?o/2) cos{aR) tan(ai?o/4)) , (12) 
2 p 

oo oo 



lim tfc = 0, / tk = +00 and > tl < og. (13) 

fe— s>oo 

k=l k=l 

then the sequence {xk)k>i generated by the algorithm satisfies 



lim p{xk, c) = 0. (14) 

fc— >oo 
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Remark 1. In practice v is given and one takes any ball B{o,R) which 
contains its support. We need the condition R < Ra- One should take R 
as small as possible for Rq and then 6 being not too small. The best choice 
is o = c and R = H{c) but they are not known a priori. If v has a finite 
support one can take for o a point of the support of v and for R the maximal 
distance from this point to another point of the support. It always works in 
a simply connected manifold of negative curvature since in this case a can 
be taken as small as we want. This is the case in our two main examples 
considered in Section 4, namely the hyperbolic space and the set of positive 
definite symmetric matrices with our specific choice of metric. Note that in 
this situation Rq and 6 can also be taken as large as we want. 

Proof: 



First we prove that for all r G [Rq, R], if Xk € B{c, r) then Xk+i € B{c, r): 
if p{xk, c) < -Ro/2 it is clear since 5 < Rq/2. If p(xfc, c) > Ro/2 we prove that 
p(xfc+i,c) < p(xfc,c). Let yk+i = 7/i-(^^)(f (x^, z/)): yk+i e supp(i/) is such 
that H{xk) = p{xk,yk+i); consider the triangle cxkyk+i- Let a = p{xk,yk+i), 
b = p{yk+i,c) and r = p{c,Xk), Xk the angle corresponding to the point Xk- 
By Alexandrov comparison theorem (in fact Corollary [T] in the appendix) Xk 
is smaller than the same in constant curvature a^. This together with the 
law of cosines in spherical geometry (Theorem 111 in the appendix) yields 



cos ah — cos ar cos aa 

COSXk > ■ 

sm ar sm aa 

Now r > Ro/2., b < r* and a > r* so 

cosar*(l — cos(Q;i?o/2)) ^ / ^ /x\ n / 7-> /.\ 

cosxfc > = cosar tan(aito/4) > cos ait tan(Q;ito/4). 

sin(aito/2) 

(15) 

Consider now the triangle cxkXk+i and let / = p(c, Xk+i)- Recall p{xk, Xk+i) = 
tk+i- Now by Toponogov theorem (Theorem [s] in the appendix) / is smaller 
than the same in constant curvature — /3^. This together with first law of 
cosines in hyperbolic geometry (Theorem 111 in the appendix) yields 



cosh f3f < cosh f3r cosh f3tk+i — cos Xk sinh f3r sinh f3tk+i (16) 



which implies by Eq. 15 



cosh/3/ < cosh(/3r) cosh/3tfc+i — cos ai?tan(ai?o/4) sinh(/5r) sinh/3tfc+i. 

(17) 
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Let us check that the condition on S imphes that the right hand side is smaller 
than cosh /3r: we want to prove 

cosh (/3r) (cosh — 1) < cosai?tan(ai?o/4) sinh(/3r) sinh/3tfc+i 

or equivalently 

cosh (3tk+i — 1 
sinh ptk+i 

But 



< cosai?tan(ai?o/4)tanh(/3r). (18) 



cosh /3tk+i - 1 _ ^^^^ f (itk+i \ 
sinh/3tfc+i V 2 y 



and tfc+i < ^7 > -Ro/2, so that Eq. 18 is implied by 

tanh — cos a;-Rtan(ai?o/4) tanh . (19) 

Now clearly the condition on 5 implies Eq. [T9j 
So we have proved that p(c, Xk+i) < p{c, Xk). 

Then we prove that there exists r] > such that if Xk G B{c, R)\B{c, Rq) 
then 

cosh (/3p(c, Xfc+i)) 



cosh (/3p(c, Xk)) 

From Eq. IlTl we obtain 



<l-Vtk+i. (20) 



< cosh f3tk+i — cos aR tan(ai?o / 4) tanh(/3r) sinh (3tk+i 

< cosh f3tk+i — cos aR tan(Q;-Ro / 4) tanh(/?i?o) sinh /3tk+i 

< 1 - 2(cosai?tan(ai?o/4)tanh(^i?o)cosh(/3tfc+i/2) 
-sinh(/3tfc+i/2))sinh(/3tfc+i/2) 

< 1 - (cosai?tan(ai?o/4) tanh(^i?o) cosh(/3tfc+i/2) - sinh(/3tfc+i/2)) jStk+i 

< 1 — (cosa-Rtan(Q;-Ro/4) tanh(^i?o) 

— cosQ;-Rtan(a_Ro/4) tanh(/3_Ro/2)) cosh.{(3tk+i/2)f3tk+i 
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where we used Eq. [12] in the last inequahty. So 
cosh /3p(c, Xfc+i) 



< 1 — (cos Q;i?tan(a/2o/4) tanh(/3/2| 







cosh /3p(c, Xk 

— cos ai?tan(ai?o/4) tanh(/3i?o/2))/?tA:+i 

(21) 

and this gives Eq. [20| 

oo 

At this stage, since ^^t/c = oo, we can conclude that there exists ko such 
fc=i 

that cosh (/3p(c, Xfcj,)) < cosh(/3_Ro) so Xko G B{c,Ro). Moreover from the 
first part of the proof we have that for all k > ko, Xk € B{c,Rq). 

Now we use the fact that on B{c,Rq), H is convex and satisfies Eq. [5| 



By boundedness of the Hessian of square distance to c (see |Yang| ( |2010[ ) 
Lemma 1.1 for details), we have for k > ko 



p^(c,Xfc+l) < 

(R -\- R \ 
— j tl+i 

(22) 

with 

C(r,/3) = 2r/3cotanh(2/3r). (23) 

Now letting yk+i = 7//(^^)(f (x^, z/)) we have H > p{-,yk+i) since yk+i E 
supp(z/). We remark that p'^{-,yk+i) is convex on B{c,Ro) by the fact that 
for all z G B{c,Rq) and y G supp(i^), p{z,y) < Ra- Moreover we have 
H{xk) = p{xk,yk+i)- As a consequence, we get 

H{c) - H{xk) > p^{c, yk+i) - p^{xk, yk+i) 
> -2 {exp~l c,jo{v{xk,iy))) 

and this implies by Lemma [T] 

- 2 (exp-^^ c, 'Joivixk, z/))) < -rp2(c, Xk). (24) 



Plugging into Eq. 22 yields 



p\c, Xk+i) < (1 - Ttk+i)p\c, Xk) + C[ ?:^1±1^ f3 ) (25) 
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We recall from here the standard argument to prove that p^(c, x^) converges 
to 0. Let 

a = limsup p^(c, Xk)- 



Iterating Eq. 25 yields for £ > 1 



p'^{c,Xk+£) < ]^(1 - Ttk+j)p^{c,Xk) +Cy^J 



2 

k+j 



with C = C ( ^"+-^ , /3) . Letting £ — )■ oo and using the fact that tk^j = oo, 
which implies 



11(1 - Tt,+j] 



0, 



we get 



a < 



■k+j- 



Finally using Yl'jLi < oo we obtain that lim^^oo J2'jLi = 0, so a = 0. 



□ 



Remark 2. In Theorem^\ it looks difficult to find a larger 6. The choice is 
almost optimal to have p{c,Xk+i) < p{c,Xk) outside B{c,Rq). On the other 



hand Eq. \21j yields an explicit value for rj in Eq. \20\ and this in turn can be 
used to find an explicit rj' > such that 



p^(c, Xk+i) < (1 - ri'tk+i)p^{c, Xk), tk+i <6 A 1/r]' 

r 



For the speed of convergence, taking t^ 



k + 1 



(26) 

we proceed as in Propo- 



sition 4.10 of Yang (2010). We use the following lemma, borrowed from the 



paper of Nedic and Bertsekas (2000): 



Lemma 2. Let {uk)k>i be a sequence of nonnegative real numbers such that 



Uk+l < 1 



A 



k + 1 
12 



Uk + 



{k + iy 



where X and ^ are positive constants. Then 



k+1 

(A-l)(fe+2) ^ (A:+2)^-i J IJ A ^ 1. 



r 

Proposition 1. Choosing tk = -, , letting such that for all k > ko, 

rC ~1~ 1 

Xk e B{c,Ro), 



(A-l){fc+2) 

where X = rr (with r given in Lemma^ and ^ = r'^C ( -^"+^ ^ /3] . 
Proof: 

This is a direct consequence of lemma |2] and inequality Eq. 25, valid for 
k>kQ. □ 

Remark 3. From the estimate of rj given by Eq. [2l\ one can get an estimate 



of ko. Another possibility is to replace r by r Arj' in Eq. 25 with r]' defined 
in Eq. 26_. Then Proposition [I] is valid for all k > 1 without the condition 



Xk e B{c,Ro). 



Remark 4. The proof of Theorem 1 works for Rq defined in Eq. 1J_. It also 
works for any smaller positive value. It is better to have Rq large so that Xk 
rapidly enters the ball B{c,Ro). On the other hand when Rq is small and 
Xk is already in this ball then one can take t close to ^cota.n{aRa) . Again 
explicit estimates are possible. 



4. Two case studies 

In order to implement algorithm GEO-ALG (a specialization of RIE-ALG 
for point clouds with step sizes = we need to describe the geodesies of 
the underlying manifold, and find an intermediate point m = Geodesic(p, q, t) 
on the geodesic passing through p and q such that p(j),m)=t p{p, q). 
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4-1- Hyperbolic manifold 

A hyperbolic manifold is a complete Riemannian ci-dimensional manifold 
of constant sectional curvature —1 that is isometric to the real hyperbolic 
space. There exists several models of hyperbolic geometry. Here, we consider 
the planar non-conformal Klein model where geodesies are straight lines. 



See Nielsen and Nock (2010). Although there exists no known closed-form 



formula for the hyperbolic centroid (p = 2), Welzl's minimax algorithm gen- 



eralizes to the Klein disk as described in Nielsen and Nock (2010) to compute 



exactly the hyperbolic 1-center. The Klein Riemannian distance on the unit 
disk is defined by 



p{p, q) = arccosh 



1 



y^{l -p^p)(l - q^q) 



(27) 



where arccosh(x) = log(x + a/x^ — 1), and the geodesic passing through p 
and q is the straight line segment 



lt{p,q) = 
Finding m such that p{p, m) 



l-t)p + tq, t G [0,1]. (28) 
=t p{p, q) cannot be solved in closed-form so- 



lution (except for t = |, see [Nielsen and Nock (2010)), so that we rather pro- 
ceed by a bisection search algorithm on parameter t up to machine precision. 
Figure [l] shows the snapshots of our implementation in Java Processing]^ 

Figure [2] plots the convergence rate of the GEO-ALG algorithm. The 
code is publicly available on-line for reproducible research. 

4-2. Manifold of symmetric positive definite matrices 

A d X d matrix M with real entries is said symmetric positive definite 
(SPD) iff. it is symmetric (M = M^), and that for all x ^ 0, x^Mx > 0. 
The set oi dx d SPD matrices forms a smooth manifold of dimension ^^i^t^_ 



We refer to Lang (1999) (Chapter 12) for a description of the geometry of 



SPD matrices. See also the work of Ji (2007) for optimization on matrix 



manifolds. The geodesic linking (matrix) point P to point Q is given by 



7t(p,g) = p^ (^p-^^Qp-^^y p"^.^ 



(29) 



processing.org 



14 






Fourth iteration 



after 104 iterations 



Figure 1: Snapshots of the GEO-ALG algorithm implemented for the hyperbolic Klein 
disk: The large black disk and the white disk denote the current center and farthest 
point, respectively. The linked path shows the trajectory of the centers as the number of 
iterations increase. On-line demo available at http : //www . inf ormationgeometry . org/] 
|RieincLnnMinim8Lx/| 
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Figure 2: Convergence rate of the GEO-ALG algorithm for the hyperbohc disk for the 
first 200 iterations. The horizontal axis denotes the number of iterations and the vertical 
axis (a) the relative Klein distance between the current center and the optimal 1-center 
(approximated for a large number of iterations), (b) the radius of the smallest enclosing 
ball anchored at the current center. 
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where the matrix function h{M) is computed from the singular value decom- 
position M = UDV^ (with U and V unitary matrices and D = diag(Ai, A^) 
a diagonal matrix of eigenvalues) as h{M) = Udia.g{h{Xi), h{\d))V^ . 
For example, the square root function of a matrix is computed as Ms = 

Udiag{^/X~l,...,^/X~d)V^. 

In this case, finding t such that 

\\\og{p-'QY\\l = r\\\ogP-'Q\\l, (30) 
where || ■ \\f denotes the Frobenius norm yields to t = r. Indeed, consider 



Ai, Xd the eigenvalues of P ^Q, then Eq. 30 amounts to find 



X: log^ AH E log^ A, = log' K. (31) 

i=l i=l i=l 

That is t = r. 

Figure |3] displays the plots of the convergence rate of the algorithm for 
the SPD manifold. 



5. Concluding remarks and discussion 



We described a generalization of the 1-center algorithm of Badoiu and 



Clarkson (2003) to arbitrary Riemannian geometry, and proved the conver- 



gence under mild assumptions. This proves the existence of Riemannian core- 
sets for optimization. This 1-center building block can be used for /c-center 
clustering. Furthermore, the algorithm can be straightforwardly extended to 
sets of geodesic balls. 

An open-source source code implementation in Java™ for reproducible 
research is available on-line at 



http : //www . inf ormationgeometry . org/RiemannMinimax/] 
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Riemannian distance between current SPD center and minimax SPD centei Radius of the smallest enclosing Riemannian ball anchored at current SPD ce 
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(a) (b) 

Figure 3: Convergence rate of the GEO-ALG algorithm for the SPD Riemannian manifold 
(dimension 5) for the first 200 iterations. The horizontal axis denotes the number of 
iterations i and the vertical axis (a) the relative Riemannian distance between the current 
center and the optimal 1-center c* , where and r* are approximated for a 

large number of iterations) , (b) the radius of the smallest enclosing SPD ball anchored 
at the current center. 
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6. Appendix: Some notions of Riemannian geometry 



In this section, we recall some basic notions of Riemannian geometry used 



throughout the paper. For a complete presentation, we refer to Cheeger and 



EbinJ(|1975j). 

We let M be a Riemannian manifold and (■,■) the Riemannian metric, 
which is a definite positive bilinear form on each tangent space T^M, and 
depends smoothly on x. The associated norm in T^M will be denoted by 
II ■ ||: ||m|| = {u,uy^'^. We denote by p{x,y) the distance between two points 
on the manifold M: 

p{x,y) = mf{[ WmUt, y^eC\[0,l],M), <^(0) = a;, ^(1) = y 



A geodesic in M is a smooth path which locally minimizes the distance 
between two points. In general such a curve does not minimize it globally. 
However it is true in all the sets we are considering in this paper. Given a 
vector V G TM with base point x, there is a unique geodesic started at x with 
speed V at time 0. It is denoted by 1 1— )■ exp^(tv) or compactly by t i— )■ 7i(w). 
It depends smoothly on v but it has in general finite lifetime. A geodesic 
defined on a time interval [a, b] is said to be minimal if it minimizes the 
distance from the image of a to the image of b. If the manifold is complete, 
taking x,y & M, there exists a minimal geodesic from a; to y in time 1. In all 
the scenarii we are considering in this paper, the minimal geodesic is unique 
and depends smoothly on x and y, and we denote it by 7.(x,?/) : [0, 1] — )■ M, 
t I— )■ 7t(a;, I/) with the conditions 7o(x, j/) = x and 7i(x,y) = y. A subset U 
of M is said to be convex if for any x,y & U, there exists a unique minimal 
geodesic 7.(x,?/) in M from x to y, this geodesic fully lies in U and depends 
smoothly on x, y, t. 

The injectivity radius of M, denoted by inj(M), is the largest r > such 
that for all x G M, the map exp^. restricted to the open ball in T^M centered 
at with radius r is an embedding. 

Given x G M, u, v two non collinear vectors in T^M, the sectional curva- 
ture Sect(M, v) = K is a number which gives information on how the geodesies 
issued from x behave near x. More precisely the image by exp^ of the circle 
centered at of radius r > in Span(M, v) has length 

27rS'x(r) + o(r^) as r 
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with 



For instance, if > 0, exp^(Span(M, is near x approximatively a 2- 

dimensional sphere with radius , In fact, if M is simply connected and 

all the sectional curvatures are equal to the same K > 0, then M is a 

dimensional sphere with radius , , where d is the dimension of M. If M 

is simply connected and all the sectional curvatures are equal to the same 
K < 0, we say that M is a d-dimensional hyperbolic space with curvature K. 

An upper bound (resp. lower bound) of sectional curvatures is a number 
a such that for all non collinear u, v in the same tangent space, Sect(ti, v) < a 
(resp. Sect(M,v) > a). In the paper, we used a positive upper bound and 
a negative lower bound — a,(3 > 0. 

The existence of the upper bound for sectional curvatures makes possi- 
ble to compare geodesic triangles, by Alexandrov theorem (see |Chavel|(l2003|)). 



Theorem 2. Let Xi, X2, X3 G M satisfy Xi 7^ X2, Xi ^ X3 and 
p{xi, X2) + p{x2, X3) + p{x3, xi) < 2 min jinjM, ^ | 

where a > is such that o? is an upper bound of sectional curvatures. Let the 
minimizing geodesic from xi to X2 and the minimizing geodesic from xi to X3 
make an angle 9 at Xi. Denoting by the 2-dimensional sphere of constant 
curvature cP' (hence of radius l/aj and p the distance in , we consider 
points xi,X2,X3 G 5*^2 such that p{xi,X2) = p{xi,X2), pixi^xs) = p{xi,xs). 
Assume that the minimizing geodesic from xi to X2 and the minimizing 
geodesic from xi to x^ also make an angle 6 at xi. 
Then we have p{x2,X3) > p{x2,Xs). 

Instead of prescribing the angle in the comparison triangle in the sphere, it 
is possible to prescribe the third distance: 

Corollary 1. The assumption are the same as in Theorem except that 
we assume that p{x2,xz) = p{x2,X3) (all the distances are equal), but the 
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minimizing geodesic from Xi to X2 and the minimizing geodesic from Xi to 
X3 now make an angle 6 at xi. 
Then we have 6 > 6. 

There also exists a comparison result in the other direction, called To- 
pogonov's theorem. 

Theorem 3. Assume (3 > is such that — is a lower bound for sectional 
curvatures in M. Let xi,X2,X3 G M satisfy Xi 7^ X2, Xi 7^ X3. Let the min- 
imizing geodesic from xi to X2 and the minimizing geodesic from xi to X3 
make an angle 9 at xi. Denoting by -^^^2 the hyperbolic 2-dimensional space 
of constant curvature — and p the distance in -f^^^2, we consider points 
xi,X2,xs G such that p{xi,X2) = p{xi,X2), p(xi,X3) = p{xi,X3). As- 

sume that the minimizing geodesic from xi to X2 and the minimizing geodesic 
from Xi to X3 also make an angle 9 at Xi. 
Then we have p{x2,Xz) < p{x2,X3). 

Triangles in the sphere and in the hyperbolic space H^^2 have explicit 
relations between distance and angles as we will see below. This combined 
with Theorems |2] and |3] and Corollary [T] allow to find related bounds in M, 
which are intensively used in our proofs. 

In this paper, we only use the first law of cosines in S^2 and in -ff^^2 (see 



e.g., the paper of Ratcliffe (1994) Theorem 2.5.3 and Theorem 3.5.3). 



Theorem 4. If 9i, 92, 9^ are the angles of a triangle in S^2 and xi, X2, X3 are 
the lengths of the opposite sides, then 

cos(q;x3) — cos(axi) cos(aa;2) 



COS^^3 



If 9i, 92, 93 are the angles of a triangle in H'^r,2 and xi,X2,X3 are the lengths 



sml^axij sm[ax2) 



of the opposite sides, then 

^ cosh(/3a;i) cosh(/3a;2) — cosh(/3x3) 
^ sinh(/3xi) sinh(/3a;2) 
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